23 February 2023
Implementing the functions for reading in genos from the rds files
created by microhaplot.
I used R-main/01-compile-and-tidy-rds-files.R to
generate an rds file with the genotype info in
data/processed.
The R function uses a minimum read depth of 10 reads for the first
allele and 6 reads for the second allele and a minimum allele balance of
0.4 (this can be modified in
R/microhaplot-genos-funcs.R).
The VCF file used to generate these rds files is
BFAL_reference_filtered.vcf, currently on Sedna. I might
return to this by merging additional samples into that VCF file, but
those would be bycatch - not reference samples - and the benefit to
doing so is uncertain.
library(tidyverse)
library(readxl)
library(stringr)
library(lubridate)
# read in rds file with genotypes
genos_long <- read_rds("../data/processed/called_genos.rds")
head(genos_long)
I need to vet both the loci and the individuals for missing data,
etc.
Locus evaluation
How many loci and how many alleles?
# alleles
genos_long %>%
filter(!is.na(allele)) %>%
select(locus, allele) %>%
unique()
# loci
genos_long %>%
filter(!is.na(allele)) %>%
select(locus, allele) %>%
unique() %>%
group_by(locus) %>%
tally() %>%
arrange(desc(n))
NA
NA
484 alleles across 190 loci with between 1-8 alleles per locus.
Missing data:
# missing data across loci
locs_to_toss <- genos_long %>%
group_by(locus) %>%
mutate(missingness = ifelse(is.na(allele), 1, 0)) %>%
summarise(sum(missingness)) %>% # 190 loci x2 = total
filter(`sum(missingness)`>190) %>% # more than 50% missing data
select(locus) # drop those loci for now and see how the assignment goes
# just the keepers
genos_locs_filtered <- genos_long %>%
anti_join(., locs_to_toss)
Joining, by = "locus"
# summary of remaining loci
genos_locs_filtered %>%
filter(!is.na(allele)) %>%
select(locus, allele) %>%
unique()
genos_locs_filtered %>%
filter(!is.na(allele)) %>%
select(locus, allele) %>%
unique() %>%
group_by(locus) %>%
tally() %>%
arrange(desc(n))
NA
398 alleles in 159 loci with 1-8 alleles per locus.
COME BACK TO THIS…
First, look at the loci for missingness and >2 haplotypes in an
individual [This might be masked by the function that reads in the rds
file??]
genos_long %>%
group_by(gtseq_run, id, locus) %>% # there should be no more than 2 alleles for a given indiv/locus
tally() %>%
filter(n > 2)
NA
Missing data in individuals
Total number of loci = 159 Total number of alleles = 318
Total number of samples = 1,055
inds_to_toss <- genos_locs_filtered %>%
group_by(gtseq_run, id) %>%
mutate(missingness = ifelse(is.na(allele), 1, 0)) %>%
summarise(sum(missingness)) %>%
arrange(desc(`sum(missingness)`)) %>%
filter(`sum(missingness)` > 63) # remove samples with >20% missing data
`summarise()` has grouped output by 'gtseq_run'. You can override using the `.groups` argument.
# just the keepers
genos_locs_ind_filtered <- genos_locs_filtered %>%
anti_join(., inds_to_toss)
Joining, by = c("gtseq_run", "id")
There’s a lot of missing data, unfortunately. We might use a higher
threshold for missing data for loci to retain more individuals if
possible.
50% missing data threshold = 993 inds to keep.
974/1055
[1] 0.9232227
Take a look at that dataset
genos_locs_ind_filtered %>%
unite(gtseq_run, id, col = "sample", remove = F) %>%
ggplot(aes(x = reorder(locus, depth), y = reorder(sample, depth), fill = log10(depth))) +
geom_tile()

self-assignment
Doing a sanity check with the reference baseline
# first make integers of the alleles
alle_idxs <- genos_locs_ind_filtered %>%
dplyr::select(gtseq_run, id, locus, gene_copy, allele) %>%
group_by(locus) %>%
mutate(alleidx = as.integer(factor(allele, levels = unique(allele)))) %>%
ungroup() %>%
arrange(gtseq_run, id, locus, alleidx) # rubias can handle NA's, so no need to change them to 0's
Add population information from metadata:
# metadata for reference samples
meta1 <- readxl::read_xlsx("../data/BFAL_WGS_01_finalPlateMap.xlsx", sheet = "combinedPlate1And2SamplePops")
samplelist <- read_csv("../data/gtseq5_samplelist.csv") %>%
mutate(gtseq_run = "gtseq5") %>%
rename(id = sample)
-- Column specification --------------------------------------------------------------------------------------------------------------------
cols(
Sample_ID = col_character(),
Sample_Plate = col_character(),
Sample_Well = col_character(),
population = col_character(),
sample = col_character()
)
metadata <- samplelist %>%
left_join(., meta1, by = c("Sample_ID" = "ind_ID")) %>%
mutate(population = ifelse(!is.na(pop), pop, population)) %>%
mutate(population = ifelse(is.na(population), "bycatch", population)) %>%
select(-pop)
# just reference samples
ref_samples <- metadata %>%
filter(!population %in% c("bycatch","NTC")) %>%
select(Sample_ID, gtseq_run, id, population) %>%
left_join(., genos_locs_ind_filtered, by = c("gtseq_run", "id"))
ref_pops <- ref_samples %>%
select(gtseq_run, id, Sample_ID, population) %>%
unique()
alle_idxs %>%
filter(gtseq_run == "gtseq5") %>%
select(id) %>%
unique() %>%
inner_join(., ref_pops) %>%
left_join(., alle_idxs)
Joining, by = "id"Joining, by = c("id", "gtseq_run")

Looks like below 90% would be a reasonable cut-off?
10 mis-assignments, but only 4 with a scaled-likelihood > 0.9 and
all of those are within the Hawaiian colonies.
Pretty good!
So overall:
129 samples in the reference baseline (10 dropped out bec of missing
data).
90% of samples correctly assigned at 90% scaled likelihood.
117/129
[1] 0.9069767
If we only look at samples that were assigned at > 90%
threshold:
117/121
[1] 0.9669421
Which would be 97% accurate assignment.
Quick look at z-scores:

Woah. There’s an outlier.
Maybe a different species accidentally? I can double-check the sample
number with metadata that Jessie sent. Other than that, things are
looking good to proceed with mixture assignment.
Mixture bycatch assignment
tmp <- bind_rows(ref_two_col, mix_two_col)
mix <- tmp %>%
filter(sample_type == "mixture")
ref <- tmp %>%
filter(sample_type == "reference")
# mixture analysis
bycatch_assign <- rubias::infer_mixture(reference = ref, mixture = mix, gen_start_col = 5)
Collating data; compiling reference allele frequencies, etc. time: 0.32 seconds
Computing reference locus specific means and variances for computing mixture z-scores time: 0.03 seconds
Working on mixture collection: bycatch with 845 individuals
calculating log-likelihoods of the mixture individuals. time: 0.08 seconds
performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.08 seconds
tidying output into a tibble. time: 0.06 seconds
Bycatch results for mixture

845 bycatch samples, most of which are assigned at high
probability.
Before going deeper into the assignment results, let’s just confirm
the loci we’re using are okay.
Output genotypes/data for evaluating the loci in
10-locus-summary-and-eval.Rmd
tmp %>%
write_csv("csv_outputs/two_column_dataset.csv")
genos_locs_ind_filtered %>%
write_rds("../data/processed/genotypes_loc_ind_filtered.rds")
LS0tDQp0aXRsZTogIjA5LXJlZmVyZW5jZS1iYXNlbGluZS1taXh0dXJlLWFuYWx5c2lzIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KMjMgRmVicnVhcnkgMjAyMw0KDQoNCg0KSW1wbGVtZW50aW5nIHRoZSBmdW5jdGlvbnMgZm9yIHJlYWRpbmcgaW4gZ2Vub3MgZnJvbSB0aGUgcmRzIGZpbGVzIGNyZWF0ZWQgYnkgbWljcm9oYXBsb3QuDQoNCkkgdXNlZCBgUi1tYWluLzAxLWNvbXBpbGUtYW5kLXRpZHktcmRzLWZpbGVzLlJgIHRvIGdlbmVyYXRlIGFuIHJkcyBmaWxlIHdpdGggdGhlIGdlbm90eXBlIGluZm8gaW4gYGRhdGEvcHJvY2Vzc2VkYC4NCg0KVGhlIFIgZnVuY3Rpb24gdXNlcyBhIG1pbmltdW0gcmVhZCBkZXB0aCBvZiAxMCByZWFkcyBmb3IgdGhlIGZpcnN0IGFsbGVsZSBhbmQgNiByZWFkcyBmb3IgdGhlIHNlY29uZCBhbGxlbGUgYW5kIGEgbWluaW11bSBhbGxlbGUgYmFsYW5jZSBvZiAwLjQgKHRoaXMgY2FuIGJlIG1vZGlmaWVkIGluIGBSL21pY3JvaGFwbG90LWdlbm9zLWZ1bmNzLlJgKS4NCg0KVGhlIFZDRiBmaWxlIHVzZWQgdG8gZ2VuZXJhdGUgdGhlc2UgcmRzIGZpbGVzIGlzIGBCRkFMX3JlZmVyZW5jZV9maWx0ZXJlZC52Y2ZgLCBjdXJyZW50bHkgb24gU2VkbmEuIEkgbWlnaHQgcmV0dXJuIHRvIHRoaXMgYnkgbWVyZ2luZyBhZGRpdGlvbmFsIHNhbXBsZXMgaW50byB0aGF0IFZDRiBmaWxlLCBidXQgdGhvc2Ugd291bGQgYmUgYnljYXRjaCAtIG5vdCByZWZlcmVuY2Ugc2FtcGxlcyAtIGFuZCB0aGUgYmVuZWZpdCB0byBkb2luZyBzbyBpcyB1bmNlcnRhaW4uDQoNCg0KDQoNCmBgYHtyIGxvYWQtbGlicmFyaWVzfQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KHJlYWR4bCkNCmxpYnJhcnkoc3RyaW5ncikNCmxpYnJhcnkobHVicmlkYXRlKQ0KYGBgDQoNCg0KYGBge3J9DQojIHJlYWQgaW4gcmRzIGZpbGUgd2l0aCBnZW5vdHlwZXMNCmdlbm9zX2xvbmcgPC0gcmVhZF9yZHMoIi4uL2RhdGEvcHJvY2Vzc2VkL2NhbGxlZF9nZW5vcy5yZHMiKQ0KDQpoZWFkKGdlbm9zX2xvbmcpDQpgYGANCg0KSSBuZWVkIHRvIHZldCBib3RoIHRoZSBsb2NpIGFuZCB0aGUgaW5kaXZpZHVhbHMgZm9yIG1pc3NpbmcgZGF0YSwgZXRjLg0KDQojIyBMb2N1cyBldmFsdWF0aW9uDQoNCkhvdyBtYW55IGxvY2kgYW5kIGhvdyBtYW55IGFsbGVsZXM/DQpgYGB7cn0NCiMgYWxsZWxlcw0KZ2Vub3NfbG9uZyAlPiUNCiAgZmlsdGVyKCFpcy5uYShhbGxlbGUpKSAlPiUNCiAgc2VsZWN0KGxvY3VzLCBhbGxlbGUpICU+JQ0KICB1bmlxdWUoKQ0KDQojIGxvY2kNCmdlbm9zX2xvbmcgJT4lDQogIGZpbHRlcighaXMubmEoYWxsZWxlKSkgJT4lDQogIHNlbGVjdChsb2N1cywgYWxsZWxlKSAlPiUNCiAgdW5pcXVlKCkgJT4lDQogIGdyb3VwX2J5KGxvY3VzKSAlPiUNCiAgdGFsbHkoKSAlPiUNCiAgYXJyYW5nZShkZXNjKG4pKQ0KDQogIA0KYGBgDQo0ODQgYWxsZWxlcyBhY3Jvc3MgMTkwIGxvY2kgd2l0aCBiZXR3ZWVuIDEtOCBhbGxlbGVzIHBlciBsb2N1cy4NCg0KTWlzc2luZyBkYXRhOg0KYGBge3J9DQojIG1pc3NpbmcgZGF0YSBhY3Jvc3MgbG9jaQ0KbG9jc190b190b3NzIDwtIGdlbm9zX2xvbmcgJT4lDQogIGdyb3VwX2J5KGxvY3VzKSAlPiUNCiAgbXV0YXRlKG1pc3NpbmduZXNzID0gaWZlbHNlKGlzLm5hKGFsbGVsZSksIDEsIDApKSAlPiUNCiAgc3VtbWFyaXNlKHN1bShtaXNzaW5nbmVzcykpICU+JSAjIDE5MCBsb2NpIHgyID0gdG90YWwNCiAgZmlsdGVyKGBzdW0obWlzc2luZ25lc3MpYD4xOTApICU+JSAjIG1vcmUgdGhhbiA1MCUgbWlzc2luZyBkYXRhDQogIHNlbGVjdChsb2N1cykgIyBkcm9wIHRob3NlIGxvY2kgZm9yIG5vdyBhbmQgc2VlIGhvdyB0aGUgYXNzaWdubWVudCBnb2VzDQoNCiMganVzdCB0aGUga2VlcGVycw0KZ2Vub3NfbG9jc19maWx0ZXJlZCA8LSBnZW5vc19sb25nICU+JQ0KICBhbnRpX2pvaW4oLiwgbG9jc190b190b3NzKQ0KDQpgYGANCmBgYHtyfQ0KIyBzdW1tYXJ5IG9mIHJlbWFpbmluZyBsb2NpDQpnZW5vc19sb2NzX2ZpbHRlcmVkICU+JQ0KICBmaWx0ZXIoIWlzLm5hKGFsbGVsZSkpICU+JQ0KICBzZWxlY3QobG9jdXMsIGFsbGVsZSkgJT4lDQogIHVuaXF1ZSgpDQoNCmdlbm9zX2xvY3NfZmlsdGVyZWQgJT4lDQogIGZpbHRlcighaXMubmEoYWxsZWxlKSkgJT4lDQogIHNlbGVjdChsb2N1cywgYWxsZWxlKSAlPiUNCiAgdW5pcXVlKCkgJT4lDQogIGdyb3VwX2J5KGxvY3VzKSAlPiUNCiAgdGFsbHkoKSAlPiUNCiAgYXJyYW5nZShkZXNjKG4pKQ0KDQpgYGANCjM5OCBhbGxlbGVzIGluIDE1OSBsb2NpIHdpdGggMS04IGFsbGVsZXMgcGVyIGxvY3VzLg0KDQoNCiMjIENPTUUgQkFDSyBUTyBUSElTLi4uDQoNCkZpcnN0LCBsb29rIGF0IHRoZSBsb2NpIGZvciBtaXNzaW5nbmVzcyBhbmQgPjIgaGFwbG90eXBlcyBpbiBhbiBpbmRpdmlkdWFsIFtUaGlzIG1pZ2h0IGJlIG1hc2tlZCBieSB0aGUgZnVuY3Rpb24gdGhhdCByZWFkcyBpbiB0aGUgcmRzIGZpbGU/P10NCg0KYGBge3J9DQpnZW5vc19sb25nICU+JQ0KICBncm91cF9ieShndHNlcV9ydW4sIGlkLCBsb2N1cykgJT4lICMgdGhlcmUgc2hvdWxkIGJlIG5vIG1vcmUgdGhhbiAyIGFsbGVsZXMgZm9yIGEgZ2l2ZW4gaW5kaXYvbG9jdXMNCiAgdGFsbHkoKSAlPiUNCiAgZmlsdGVyKG4gPiAyKQ0KDQpgYGANCg0KDQojIyBNaXNzaW5nIGRhdGEgaW4gaW5kaXZpZHVhbHMNCg0KVG90YWwgbnVtYmVyIG9mIGxvY2kgPSAxNTkNClRvdGFsIG51bWJlciBvZiBhbGxlbGVzID0gMzE4DQoNClRvdGFsIG51bWJlciBvZiBzYW1wbGVzID0gMSwwNTUNCg0KYGBge3J9DQppbmRzX3RvX3Rvc3MgPC0gZ2Vub3NfbG9jc19maWx0ZXJlZCAlPiUNCiAgZ3JvdXBfYnkoZ3RzZXFfcnVuLCBpZCkgJT4lDQogIG11dGF0ZShtaXNzaW5nbmVzcyA9IGlmZWxzZShpcy5uYShhbGxlbGUpLCAxLCAwKSkgJT4lDQogIHN1bW1hcmlzZShzdW0obWlzc2luZ25lc3MpKSAlPiUNCiAgYXJyYW5nZShkZXNjKGBzdW0obWlzc2luZ25lc3MpYCkpICU+JQ0KICBmaWx0ZXIoYHN1bShtaXNzaW5nbmVzcylgID4gNjMpICMgcmVtb3ZlIHNhbXBsZXMgd2l0aCA+MjAlIG1pc3NpbmcgZGF0YQ0KDQojIGp1c3QgdGhlIGtlZXBlcnMNCmdlbm9zX2xvY3NfaW5kX2ZpbHRlcmVkIDwtIGdlbm9zX2xvY3NfZmlsdGVyZWQgJT4lDQogIGFudGlfam9pbiguLCBpbmRzX3RvX3Rvc3MpDQoNCmBgYA0KVGhlcmUncyBhIGxvdCBvZiBtaXNzaW5nIGRhdGEsIHVuZm9ydHVuYXRlbHkuIFdlIG1pZ2h0IHVzZSBhIGhpZ2hlciB0aHJlc2hvbGQgZm9yIG1pc3NpbmcgZGF0YSBmb3IgbG9jaSB0byByZXRhaW4gbW9yZSBpbmRpdmlkdWFscyBpZiBwb3NzaWJsZS4NCg0KNTAlIG1pc3NpbmcgZGF0YSB0aHJlc2hvbGQgPSA5OTMgaW5kcyB0byBrZWVwLg0KDQpgYGB7cn0NCjk3NC8xMDU1DQpgYGANClRha2UgYSBsb29rIGF0IHRoYXQgZGF0YXNldA0KYGBge3J9DQpnZW5vc19sb2NzX2luZF9maWx0ZXJlZCAgJT4lDQogIHVuaXRlKGd0c2VxX3J1biwgaWQsIGNvbCA9ICJzYW1wbGUiLCByZW1vdmUgPSBGKSAlPiUNCiAgZ2dwbG90KGFlcyh4ID0gcmVvcmRlcihsb2N1cywgZGVwdGgpLCB5ID0gcmVvcmRlcihzYW1wbGUsIGRlcHRoKSwgZmlsbCA9IGxvZzEwKGRlcHRoKSkpICsNCiAgZ2VvbV90aWxlKCkNCg0KYGBgDQoNCiMjIHNlbGYtYXNzaWdubWVudCANCg0KRG9pbmcgYSBzYW5pdHkgY2hlY2sgd2l0aCB0aGUgcmVmZXJlbmNlIGJhc2VsaW5lDQoNCmBgYHtyfQ0KIyBmaXJzdCBtYWtlIGludGVnZXJzIG9mIHRoZSBhbGxlbGVzDQphbGxlX2lkeHMgPC0gZ2Vub3NfbG9jc19pbmRfZmlsdGVyZWQgJT4lIA0KICBkcGx5cjo6c2VsZWN0KGd0c2VxX3J1biwgaWQsIGxvY3VzLCBnZW5lX2NvcHksIGFsbGVsZSkgJT4lDQogIGdyb3VwX2J5KGxvY3VzKSAlPiUNCiAgbXV0YXRlKGFsbGVpZHggPSBhcy5pbnRlZ2VyKGZhY3RvcihhbGxlbGUsIGxldmVscyA9IHVuaXF1ZShhbGxlbGUpKSkpICU+JQ0KICB1bmdyb3VwKCkgJT4lDQogIGFycmFuZ2UoZ3RzZXFfcnVuLCBpZCwgbG9jdXMsIGFsbGVpZHgpICMgcnViaWFzIGNhbiBoYW5kbGUgTkEncywgc28gbm8gbmVlZCB0byBjaGFuZ2UgdGhlbSB0byAwJ3MNCg0KYGBgDQoNCg0KQWRkIHBvcHVsYXRpb24gaW5mb3JtYXRpb24gZnJvbSBtZXRhZGF0YToNCmBgYHtyfQ0KIyBtZXRhZGF0YSBmb3IgcmVmZXJlbmNlIHNhbXBsZXMNCm1ldGExIDwtIHJlYWR4bDo6cmVhZF94bHN4KCIuLi9kYXRhL0JGQUxfV0dTXzAxX2ZpbmFsUGxhdGVNYXAueGxzeCIsIHNoZWV0ID0gImNvbWJpbmVkUGxhdGUxQW5kMlNhbXBsZVBvcHMiKQ0KDQpzYW1wbGVsaXN0IDwtIHJlYWRfY3N2KCIuLi9kYXRhL2d0c2VxNV9zYW1wbGVsaXN0LmNzdiIpICU+JQ0KICBtdXRhdGUoZ3RzZXFfcnVuID0gImd0c2VxNSIpICU+JQ0KICByZW5hbWUoaWQgPSBzYW1wbGUpDQoNCm1ldGFkYXRhIDwtIHNhbXBsZWxpc3QgJT4lDQogIGxlZnRfam9pbiguLCBtZXRhMSwgYnkgPSBjKCJTYW1wbGVfSUQiID0gImluZF9JRCIpKSAlPiUNCiAgbXV0YXRlKHBvcHVsYXRpb24gPSBpZmVsc2UoIWlzLm5hKHBvcCksIHBvcCwgcG9wdWxhdGlvbikpICU+JQ0KICBtdXRhdGUocG9wdWxhdGlvbiA9IGlmZWxzZShpcy5uYShwb3B1bGF0aW9uKSwgImJ5Y2F0Y2giLCBwb3B1bGF0aW9uKSkgJT4lDQogIHNlbGVjdCgtcG9wKQ0KDQojIGp1c3QgcmVmZXJlbmNlIHNhbXBsZXMNCnJlZl9zYW1wbGVzIDwtIG1ldGFkYXRhICU+JQ0KICBmaWx0ZXIoIXBvcHVsYXRpb24gJWluJSBjKCJieWNhdGNoIiwiTlRDIikpICU+JQ0KICBzZWxlY3QoU2FtcGxlX0lELCBndHNlcV9ydW4sIGlkLCBwb3B1bGF0aW9uKSAlPiUNCiAgbGVmdF9qb2luKC4sIGdlbm9zX2xvY3NfaW5kX2ZpbHRlcmVkLCBieSA9IGMoImd0c2VxX3J1biIsICJpZCIpKQ0KICANCnJlZl9wb3BzIDwtIHJlZl9zYW1wbGVzICU+JQ0KICBzZWxlY3QoZ3RzZXFfcnVuLCBpZCwgU2FtcGxlX0lELCBwb3B1bGF0aW9uKSAlPiUNCiAgdW5pcXVlKCkNCmBgYA0KDQpgYGB7cn0NCmFsbGVfaWR4cyAlPiUNCiAgZmlsdGVyKGd0c2VxX3J1biA9PSAiZ3RzZXE1IikgJT4lDQogIHNlbGVjdChpZCkgJT4lDQogIHVuaXF1ZSgpICU+JQ0KICBpbm5lcl9qb2luKC4sIHJlZl9wb3BzKSAlPiUNCiAgbGVmdF9qb2luKC4sIGFsbGVfaWR4cykNCg0KYGBgDQoNCg0KDQpgYGB7cn0NCiMgZm9ybWF0IGZvciBydWJpYXMNCnJlZmVyZW5jZSA8LSBhbGxlX2lkeHMgJT4lDQogIGlubmVyX2pvaW4oLiwgcmVmX3BvcHMpICU+JQ0KICBzZWxlY3QoLWFsbGVsZSwgLWd0c2VxX3J1biwgLWlkKSAlPiUNCiAgc2VsZWN0KHBvcHVsYXRpb24sIFNhbXBsZV9JRCwgZXZlcnl0aGluZygpKSAlPiUNCiAgcmVuYW1lKGNvbGxlY3Rpb24gPSBwb3B1bGF0aW9uLCBpbmRpdiA9IFNhbXBsZV9JRCkNCg0KIyBtYWtlIHR3by1jb2wgZm9ybWF0DQpyZWZfdHdvX2NvbCA8LSByZWZlcmVuY2UgJT4lDQogIHVuaXRlKCJsb2MiLCAzOjQsIHNlcCA9ICIuIikgJT4lDQogIHBpdm90X3dpZGVyKG5hbWVzX2Zyb20gPSBsb2MsIHZhbHVlc19mcm9tID0gYWxsZWlkeCkgJT4lDQogIG11dGF0ZShyZXB1bml0ID0gY29sbGVjdGlvbikgJT4lDQogIG11dGF0ZShzYW1wbGVfdHlwZSA9ICJyZWZlcmVuY2UiKSAlPiUNCiAgc2VsZWN0KHNhbXBsZV90eXBlLCByZXB1bml0LCBjb2xsZWN0aW9uLCBldmVyeXRoaW5nKCkpDQoNCmBgYA0KDQoNCmBgYHtyfQ0KIyBzZWxmLWFzc2lnbm1lbnQNCmJhc2VsaW5lX2Fzc2lnbiA8LSBydWJpYXM6OnNlbGZfYXNzaWduKHJlZl90d29fY29sLCBnZW5fc3RhcnRfY29sID0gNSkNCg0KYGBgDQoNCmBgYHtyfQ0KdG9wX2Fzc2lnbiA8LSBiYXNlbGluZV9hc3NpZ24gJT4lDQogIGdyb3VwX2J5KGluZGl2KSAlPiUNCiAgc2xpY2VfbWF4KC4sIG9yZGVyX2J5ID0gc2NhbGVkX2xpa2VsaWhvb2QpIA0KDQp0b3BfYXNzaWduICU+JSAjIHRvcCBhc3NpZ25tZW50IGZvciBlYWNoIHNhbXBsZQ0KICBnZ3Bsb3QoYWVzKHggPSBzY2FsZWRfbGlrZWxpaG9vZCkpICsNCiAgZ2VvbV9oaXN0b2dyYW0oKQ0KDQpgYGANCkxvb2tzIGxpa2UgYmVsb3cgOTAlIHdvdWxkIGJlIGEgcmVhc29uYWJsZSBjdXQtb2ZmPw0KDQpgYGB7cn0NCnRvcF9hc3NpZ24gJT4lDQogIGZpbHRlcihyZXB1bml0ICE9IGluZmVycmVkX2NvbGxlY3Rpb24sIA0KICAgICAgICAgc2NhbGVkX2xpa2VsaWhvb2QgPiAwLjkpICU+JQ0KICBnZ3Bsb3QoYWVzKHggPSBzY2FsZWRfbGlrZWxpaG9vZCkpICsNCiAgZ2VvbV9oaXN0b2dyYW0oKQ0KDQpgYGANCjEwIG1pcy1hc3NpZ25tZW50cywgYnV0IG9ubHkgNCB3aXRoIGEgc2NhbGVkLWxpa2VsaWhvb2QgPiAwLjkgYW5kIGFsbCBvZiB0aG9zZSBhcmUgd2l0aGluIHRoZSBIYXdhaWlhbiBjb2xvbmllcy4NCg0KUHJldHR5IGdvb2QhDQoNClNvIG92ZXJhbGw6DQoNCmBgYHtyfQ0KdG9wX2Fzc2lnbiAlPiUNCiAgZmlsdGVyKHNjYWxlZF9saWtlbGlob29kID4gMC45ICYNCiAgICAgICAgIGNvbGxlY3Rpb24gPT0gaW5mZXJyZWRfY29sbGVjdGlvbikNCmBgYA0KMTI5IHNhbXBsZXMgaW4gdGhlIHJlZmVyZW5jZSBiYXNlbGluZSAoMTAgZHJvcHBlZCBvdXQgYmVjIG9mIG1pc3NpbmcgZGF0YSkuDQoNCjkwJSBvZiBzYW1wbGVzIGNvcnJlY3RseSBhc3NpZ25lZCBhdCA5MCUgc2NhbGVkIGxpa2VsaWhvb2QuDQpgYGB7cn0NCjExNy8xMjkNCmBgYA0KSWYgd2Ugb25seSBsb29rIGF0IHNhbXBsZXMgdGhhdCB3ZXJlIGFzc2lnbmVkIGF0ID4gOTAlIHRocmVzaG9sZDoNCmBgYHtyfQ0KdG9wX2Fzc2lnbiAlPiUNCiAgZmlsdGVyKHNjYWxlZF9saWtlbGlob29kID4gMC45KQ0KYGBgDQpgYGB7cn0NCjExNy8xMjENCmBgYA0KV2hpY2ggd291bGQgYmUgOTclIGFjY3VyYXRlIGFzc2lnbm1lbnQuDQoNClF1aWNrIGxvb2sgYXQgei1zY29yZXM6DQoNCmBgYHtyfQ0KdG9wX2Fzc2lnbiAlPiUNCiAgZ2dwbG90KGFlcyh6X3Njb3JlKSkgKw0KICBnZW9tX2hpc3RvZ3JhbSgpDQoNCg0KYGBgDQpXb2FoLiBUaGVyZSdzIGFuIG91dGxpZXIuDQoNCmBgYHtyfQ0KdG9wX2Fzc2lnbiAlPiUNCiAgZmlsdGVyKHpfc2NvcmUgPCAtMykNCg0KYGBgDQpNYXliZSBhIGRpZmZlcmVudCBzcGVjaWVzIGFjY2lkZW50YWxseT8gSSBjYW4gZG91YmxlLWNoZWNrIHRoZSBzYW1wbGUgbnVtYmVyIHdpdGggbWV0YWRhdGEgdGhhdCBKZXNzaWUgc2VudC4gT3RoZXIgdGhhbiB0aGF0LCB0aGluZ3MgYXJlIGxvb2tpbmcgZ29vZCB0byBwcm9jZWVkIHdpdGggbWl4dHVyZSBhc3NpZ25tZW50Lg0KDQoNCiMjIE1peHR1cmUgYnljYXRjaCBhc3NpZ25tZW50DQoNCg0KYGBge3J9DQojIGdldCB0aGUgZm9ybWF0IHJpZ2h0DQptaXhfdHdvX2NvbCA8LSBhbGxlX2lkeHMgJT4lDQogIGFudGlfam9pbiguLCByZWZfcG9wcykgJT4lDQogIHVuaXRlKGd0c2VxX3J1biwgaWQsIGNvbCA9ICJpbmRpdiIsIHNlcCA9ICJfIikgJT4lDQogIHNlbGVjdCgtYWxsZWxlKSAlPiUNCiAgdW5pdGUoImxvYyIsIDI6Mywgc2VwID0gIi4iKSAlPiUNCiAgcGl2b3Rfd2lkZXIobmFtZXNfZnJvbSA9IGxvYywgdmFsdWVzX2Zyb20gPSBhbGxlaWR4KSAlPiUNCiAgbXV0YXRlKHJlcHVuaXQgPSBOQSkgJT4lDQogIG11dGF0ZShzYW1wbGVfdHlwZSA9ICJtaXh0dXJlIikgJT4lDQogIG11dGF0ZShjb2xsZWN0aW9uID0gImJ5Y2F0Y2giKSAlPiUNCiAgc2VsZWN0KHNhbXBsZV90eXBlLCByZXB1bml0LCBjb2xsZWN0aW9uLCBldmVyeXRoaW5nKCkpDQoNCm1peF90d29fY29sJHJlcHVuaXQgPC0gYXMuY2hhcmFjdGVyKG1peF90d29fY29sJHJlcHVuaXQpDQpoZWFkKG1peF90d29fY29sKQ0KDQpgYGANCg0KYGBge3J9DQojIGZvciB3aGF0ZXZlciByZWFzb24sIHJ1YmlhcyBpcyBiZWluZyBmaW5pY2t5IGFib3V0IHRoZSBmb3JtYXQgYW5kIGl0J3MgZWFzaWVzdCB0byBkbyB0aGlzDQojIHRvIG1ha2Ugc3VyZSBib3RoIGRmIGFyZSBjb25zaXN0ZW50DQp0bXAgPC0gYmluZF9yb3dzKHJlZl90d29fY29sLCBtaXhfdHdvX2NvbCkNCg0KbWl4IDwtIHRtcCAlPiUNCiAgZmlsdGVyKHNhbXBsZV90eXBlID09ICJtaXh0dXJlIikNCg0KcmVmIDwtIHRtcCAlPiUNCiAgZmlsdGVyKHNhbXBsZV90eXBlID09ICJyZWZlcmVuY2UiKQ0KYGBgDQoNCg0KYGBge3J9DQojIG1peHR1cmUgYW5hbHlzaXMNCmJ5Y2F0Y2hfYXNzaWduIDwtIHJ1Ymlhczo6aW5mZXJfbWl4dHVyZShyZWZlcmVuY2UgPSByZWYsIG1peHR1cmUgPSBtaXgsIGdlbl9zdGFydF9jb2wgPSA1KQ0KDQpgYGANCg0KQnljYXRjaCByZXN1bHRzIGZvciBtaXh0dXJlDQpgYGB7cn0NCm1peF9hc3NpZ25lZCA8LSBieWNhdGNoX2Fzc2lnbiRpbmRpdl9wb3N0ZXJpb3JzICU+JQ0KICBncm91cF9ieShpbmRpdikgJT4lDQogIHNsaWNlX21heCguLCBvcmRlcl9ieSA9IFBvZlopDQoNCg0KIyBkaXN0cmlidXRpb24gb2YgUG9mWiBmb3IgdG9wIGFzc2lnbm1lbnRzDQptaXhfYXNzaWduZWQgJT4lDQogIGdncGxvdChhZXMoeCA9IFBvZlopKSArDQogIGdlb21faGlzdG9ncmFtKCkNCg0KYGBgDQo4NDUgYnljYXRjaCBzYW1wbGVzLCBtb3N0IG9mIHdoaWNoIGFyZSBhc3NpZ25lZCBhdCBoaWdoIHByb2JhYmlsaXR5Lg0KDQoNCmBgYHtyfQ0KbWl4X2Fzc2lnbmVkICU+JQ0KICBmaWx0ZXIoUG9mWiA+IDAuOSkNCg0KYGBgDQoNCkJlZm9yZSBnb2luZyBkZWVwZXIgaW50byB0aGUgYXNzaWdubWVudCByZXN1bHRzLCBsZXQncyBqdXN0IGNvbmZpcm0gdGhlIGxvY2kgd2UncmUgdXNpbmcgYXJlIG9rYXkuDQoNCg0KT3V0cHV0IGdlbm90eXBlcy9kYXRhIGZvciBldmFsdWF0aW5nIHRoZSBsb2NpIGluIGAxMC1sb2N1cy1zdW1tYXJ5LWFuZC1ldmFsLlJtZGANCg0KYGBge3J9DQojIHNhdmUgb3V0cHV0cw0KdG1wICU+JQ0KICB3cml0ZV9jc3YoImNzdl9vdXRwdXRzL3R3b19jb2x1bW5fZGF0YXNldC5jc3YiKQ0KDQpnZW5vc19sb2NzX2luZF9maWx0ZXJlZCAlPiUNCiAgd3JpdGVfcmRzKCIuLi9kYXRhL3Byb2Nlc3NlZC9nZW5vdHlwZXNfbG9jX2luZF9maWx0ZXJlZC5yZHMiKQ0KDQphbGxlX2lkeHMgJT4lDQogIHdyaXRlX3JkcygiLi4vZGF0YS9wcm9jZXNzZWQvYWxsZV9pZHhzLnJkcyIpDQoNCnJlZmVyZW5jZSAlPiUNCiAgd3JpdGVfcmRzKCIuLi9kYXRhL3Byb2Nlc3NlZC9yZWZlcmVuY2VfZ2Vub3MucmRzIikNCmBgYA0KDQo=